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ABSTRACT 

We discuss the effect of chemical separation as matter freezes at the base of the ocean of an accreting 
neutron star, and argue that the retention of light elements in the liquid acts as a source of buoyancy 
that drives a slow but continual mixing of the ocean, enriching it substantially in light elements, and 
leading to a relatively uniform composition with depth. We first consider the timescales associated with 
different processes that can redistribute elements in the ocean, including convection, sedimentation, 
crystallization, and diffusion. We then calculate the steady state structure of the ocean of a neutron 
star for an illustrative model in which the accreted hydrogen and helium burns to produce a mixture 
of O and Se. Even though the H/He burning produces only 2% oxygen by mass, the steady state ocean 
has an oxygen abundance more than ten times larger, almost 40% by mass. Furthermore, we show 
that the convective motions transport heat inwards, with a flux of ~ 0.2 MeV per nucleon for an O-Se 
ocean, heating the ocean and steepening the outwards temperature gradient. The enrichment of light 
elements and heating of the ocean due to compositionally-driven convection likely have important 
implications for carbon ignition models of superbursts. 

Subject headings: dense matter — stars: neutron — X-rays: binaries — X-rays: individual 



1. INTRODUCTION 

The ocean of an accreting neutron star is composed 
of a variety of elements with atomic number Z = 6 and 
larger, formed by nuclear burning of the accreted hydro- 
gen and helium at low densities. The term ocean refers 
to the fact that the Coulomb interaction energy between 
ions is greater than the thermal energy, such that the ions 
behave like a liquid. The ocean is of interest as the site of 
long duration thermonuclear flashes such as superbursts 
(Cumming & Bildsten 2001; Strohmayer & Brown 2002; 
Kuulkcrs 2004) and intermediate duration bursts (in 't 
Zand et al. 2005; Cumming et al. 2006), non-radial oscil- 
lations (Bildsten & Cutler 1995; Piro & Bildsten 2005), 
and because the matter in the ocean eventually solidi- 
fies as it is compressed to higher densities by ongoing 
accretion, and so determines the thermal, mechanical 
and compositional properties of the neutron star crust 
(Haensel & Zdunik 1990; Brown k Bildsten 1998; Schatz 
et al. 1999). 

At the base of the ocean, matter freezes as it is com- 
pressed by continuing accretion, becoming part of the 
solid crust. Horowitz et al. (2007) carried out molecular 
dynamics simulations of the freezing of a mixture of 17 
species taken from a calculation of rp-process hydrogen 
and helium burning and hence representative of the kind 
of mixture expected to make up the ocean of an accret- 
ing neutron star (Schatz et al. 2001; Gupta et al. 2007). 
They found that this mixture underwent chemical sep- 
aration during crystallization, such that light elements 
(charge number Z < 20) were preferentially left behind 
in the liquid, whereas heavier elements were preferen- 
tially incorporated into the solid. In a previous paper 
(Medin & Cumming 2010, hereafter Paper I) we showed 
that this result can be understood by generalizing previ- 
ous work using fits to the free energies of the liquid and 
solid states of binary and tertiary plasmas. 

In this paper, we address the implications of chemi- 



cal separation for the structure and composition of the 
ocean. Horowitz et al. (2007) raised the question of what 
the steady-state ocean would look like, since the matter 
entering the crust is enriched in certain elements com- 
pared to others, and therefore different from the mean 
ocean composition. We investigate this question here, 
and argue that the retention of light elements in the liq- 
uid acts as a source of buoyancy that drives a slow but 
continual mixing of the ocean, enriching it substantially 
in light elements and leading to a relatively uniform com- 
position with depth. The steady state arises as the ocean 
enriches in light elements to the point where the compo- 
sition of the solid that forms upon freezing matches the 
composition of matter entering the top of the ocean. 

One motivation for studying this problem comes from 
models for superbursts which involve thermally-unstable 
carbon burning in the deep ocean of the neutron star 
(Cumming & Bildsten 2001; Strohmayer & Brown 2002). 
The energy release in these very long duration thermonu- 
clear flashes, inferred from fitting their lightcurves (Cum- 
ming et al. 2006), corresponds to carbon fractions of 
« 20%. This has been challenging to produce in models 
of the nuclear burning of the accreted hydrogen and he- 
lium. If the hydrogen and helium burn unstably, the 
amount of carbon produced is < 1% (Woosley et al. 
2004), and whereas stable burning can produce large car- 
bon fractions (Schatz et al. 2003), time-dependent mod- 
els do not show stable burning at the « 10% Eddington 
accretion rates of superburst sources (although observa- 
tionally, superburst sources show evidence that much of 
the accreted material may not burn in Type I bursts; 
in 't Zand et al. 2003). 

Perhaps even more problematic than making enough 
carbon is that carbon ignition models for superbursts 
require large ocean temperatures w 6 x 10 8 K at the ig- 
nition depth, which are difficult to achieve in standard 
models of crust heating (e.g., Cumming et al. 2006; Keek 
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et al. 2008). Similarly, Brown & Cumming (2009) in- 
ferred a large inwards heat flux in the outer crust of the 
transiently-accreting neutron stars MXB 1659-29 and KS 
1731-260 by fitting their cooling curves in quiescence. 
Both of these observations imply an additional heating 
source in the outer crust or ocean is needed. In this pa- 
per, we begin to address the question of to what extent 
chemical separation could enrich the ocean in carbon and 
other light elements, or provide a heat source that could 
alleviate some of the difficulty of matching the observa- 
tions of superbursts and transient cooling. 

We begin in <j2] by reviewing the physics of chemical 
separation, and discussing the timescales on which accre- 
tion, crystallization, diffusion, sedimentation, and con- 
vection occur, leading us to a picture of compositionally- 
driven convection. In we calculate the structure of 
the steady-state ocean for two simplified models: first, 
accretion of a two-component mixture composed of Se 
and either O or Fe; and second, accretion of a mixture 
of H and He which then burns to produce these heavier- 
element mixtures. In Sj4] we consider the effect of the 
mixing on the thermal profile, and calculate the heat- 
ing of the ocean due to the convective transport of light 
elements outwards. Finally, in $5] we discuss the impli- 
cations of our results. 

2. MIXING PROCESSES IN THE OCEAN 

2.1. Phase diagrams and chemical separation 

The degree of chemical separation on freezing can be 
understood from the phase diagram for the mixture. Fig- 
ure [1] shows two examples of phase diagrams, calculated 
as described in Paper I. The upper panel is for a mixture 
of 56 Fe (Z = 26) and 79 Se (Z = 34), the lower panel 
for a mixture of 16 (Z = 8) and 79 Sc (Z = 34). In 
the first case the diagram is of azeotrope type; in the 
second case, for which the ratio of atomic numbers Z is 
greater, the diagram is a more complicated eutectic type. 
In each case, the x-axis shows the number fraction of Se 
and the y-axis shows the inverse Coulomb coupling pa- 
rameter for the light species (Fe or O respectively) 
The Coulomb coupling parameter for species i is 



? 5/3„2 



1/3 



k B T \ 3m p ) 



1/3 



(1) 



where Z% is the charge of the ion, Y e = (Z)/(A), (Z) 
and (A) are the average charge and mass per ion of 
the mixture, pg = p/(10 9 g/cm 3 ) the density, and T 8 = 
T/(10 8 K) the temperature. For a single species of ion, 
solidification occurs when Ti > T m w 175 (e.g., Potekhin 
& Chabrier 2000). Note that for a given p and T (or 
depth in the star), Ti is nearly constant with composi- 
tion. As a fluid element is compressed by accretion to 
higher density, Ti increases, moving down in the phase 
diagram. The shaded regions represent unstable regions 
of the phase diagram. A fluid element with composition 
and Fi that lies inside the unstable region will undergo 
phase separation, separating into two phases with com- 
positions on each side of the unstable region. In this way, 
chemical separation occurs. Note that the curves that 
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Fig. 1. — (Color online) The phase diagram for crystallization of 
an 56 Fe- 79 Se mixture (top panel) and an 16 0- 79 Se mixture (bottom 
panel) in a T = 3 X 10 8 K ocean. The Coulomb coupling constants 
rp o are given in terms of r cr ; t fa 175, the value at which a single- 
species plasma crystallizes. The stable liquid region of each phase 
diagram is labeled 'L', the stable solid region(s) are labeled 'S' or 
'SI' and 'S2', and the unstable region is filled with plus symbols. 
Additionally, in each panel the composition at the top of the ocean 
is marked by a vertical dashed line, the ocean-crust boundary is 
marked by a horizontal dotted line, the composition of the liquid 
at the base of the ocean is marked by a filled square, and the 
composition of the solid(s) in the outer crust are marked by filled 
circles. 



bound the unstable region are commonly referred to as 
the liquidus and solidus curves, respectively. A liquid 
with composition of the liquidus curve at a given Ti is 
in equilibrium with the solid which has the composition 
of the solidus curve at that Ti. 

Now consider a particular mixture of Fe and Se enter- 
ing the top of the ocean with xso = 0.7 (this corresponds 
to 77% Se and 23% Fe by mass), indicated by the verti- 
cal dashed line in the upper panel of Fig. [TJ In steady 
state, the solid forming at the top of the crust must have 
this same composition, so that the freezing point must 
lie at r m /TFc ~ 1-1 as marked by the filled circle in 
Fig. [T] The corresponding liquid composition in equilib- 
rium with the solid at this T-p e is indicated by the solid 
square. The phase diagram shows that the liquid at the 
base of the ocean must have a composition xg c 0.56 
(63% Se by mass) in order to make the solid composition 
demanded by steady state. 
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The lower panel of Fig. [T] shows a more complicated 
example. The vertical dotted line marks an incoming 
composition with xs c = 0.9 (corresponding to 98% Se 
and 2% O by mass). In this case, there is no single 
solid phase with this composition. Instead, at T m /TQ ss 
3.5, a mixture of two solid phases indicated by the filled 
circles forms, and the liquid at the base of the ocean 
has xsc ~ 0.25 (63% Se by mass). Again, in order to 
reach steady state, the base of the ocean must adjust 
its composition until it is significantly enriched in light 
elements compared to the composition at the top of the 



2.2. Crystallization of solid particles 

To see how the ocean is able to adjust its composition 
profile to achieve steady state, we first note that solid 
particles crystallize and sediment out rapidly compared 
to the accretion timescale on which matter is compressed. 
The accretion time is i accr = Hp / u aC cr — y/m or 
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where Hp — P/pg = y/p is the pressure scale height, 
y is the column depth, m is the local accretion rate 
per unit area and in the last step we use the fact that 
the pressure in the ocean is dominated by relativis- 
tic degenerate electrons. We scale to a typical rate 
m 4 = m/ (10 4 g cm~ 2 s _1 ) (to set the scale note that the 
Eddington rate is rnEdd ~ 10 5 g cm~ 2 s _1 ). We choose 
the gravity gn — g/(10 14 cm/s 2 ) = 2.45 corresponding 
to a 1.4 Mq, R — 10 km neutron star. 

First, consider crystallization. The nucleation rate, the 
rate at which solid clusters of a size large enough to be 
stable are formed, is not currently well understood. The- 
ory typically predicts nucleation rates orders of magni- 
tude smaller than experiment or simulations indicate (see 
Vehkamaki 2006). For an OCP for example, the nucle- 
ation rate at Ti = 300 from the theoretical model of 
Ichimaru et al. (1983) is ~ 10~ 6 of its value from the sim- 
ulation of Daligault (2006b). Cooper & Bildsten (2008) 
recently derived a value within a factor of 10 of the sim- 
ulations, but their formulation predicts the formation of 
stable solid clusters even at Ti < 175 when the liquid 
phase should be absolutely stable with respect to the 
solid. Taking the Ichimaru et al. (1983) results, we find 
that the amount of undercooling necessary for the nucle- 
ation rate to become comparable to the accretion rate is 
Ti — r m ~ 0.ir m , but given that the theory underpre- 
dicts the simulation results, the amount of undercooling 
required is probably much less than this. Additionally, 
at the base of the ocean solid clusters do not need to 
wait for nucleation sites to form but can crystallize on 
the existing crust. 

Once a cluster forms it increases in size at the crys- 



tallization velocity which is of order v c 
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Kelton, Greer, & Thompson 1983), where a is the mean 
ion spacing and ui p is the ion plasma frequency given by 



= ^/Wp(Y e e/m p ) 



1.4 x 10 19 pg /2 (Y e /0.43) rad/s. 
In a multicomponent plasma crystal growth is slower, 
since as a solid cluster grows chemical separation means 
that the liquid surrounding the cluster is depleted more 
and more of the particles necessary to form the solid. 



The crystallization rate therefore depends on the rate at 
which diffusion can replenish the depleted particles. In 

a liquid with 100 < T < 300, where V = (Z 5 ^)Ti/zf /3 
is the average Coulomb coupling constant, the diffusion 
coefficient for species i is 
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(Horowitz et al. 2010; see also Daligault & Murillo 
2005). We assume that to bind to a cluster, new ma- 
terial must travel a distance I equal to the current size 
of the cluster; in this case, the cluster growth time is 
l 2 /D « 10- 4 Wp 1 (Va) 2 ( r / 1 00) 2 5 . For a solid cluster of 

N s particles, (I /a) 2 w A s 2/3 . 

Therefore as soon as the mixture encounters the liq- 
uidus line in the phase diagram, we expect solid parti- 
cles to rapidly form and grow. Once formed, the solid 
particles will quickly sediment out. We estimate the sed- 
imentation velocity for solid clusters following Bildsten & 
Hall (2001) and Brown et al. (2002). The sedimentation 
velocity is given by the Einstein relation 
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where AY e /Y e = 1 — Y e ^jY e ^ is the contrast between 
Y e for the solid particles and the background fluid. 
Note that in Eq. (0| we have neglected the contribu- 
tion of the ions to the buoyancy force. Usually this 
amounts to a ~ 10% correction; but when AY e — 
this is the dominant term (Mochkovitch 1983). We 
use the Stokes-Einstein relation to estimate the mobil- 
ity D/(k B T) = (47T77aA s 1/3 ) , where the OCP shear 

viscosity is 77 pa 0.2/5w p a 2 (r/100) 1 ' 5 g cm -1 s _1 (Donko 
& Nyfri 2000; Daligault 2006a). We find 
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(5) 

Comparing with the accretion velocity 7> a ccr = m/p = 
10~ 5 rri4 : pg 1 cm/s, we find that the critical cluster size 
above which the particles sediment out is 

AT OOO -0.6 f T &\ • 3/2 ( AYe/Ye X 

A s cr it = 200 p Q — mJ 

' H9 V 3 J 4 V 0.01 

(6) 

Note that the simple nature of our estimates means that 
there are considerable uncertainties in the crystallization 
and sedimentation velocities we find here. Despite this, 
however, the timescale for a cluster to grow to a size 
A S!Cr it is so short that the conclusion that the solid par- 
ticles rapidly fall out of the ocean seems inescapable. 

Setting A S!C1 it = 1 in Eq. ([6]) provides an estimate 
of when relative separation of light and heavy elements 
(without forming solid clusters) is expected to occur. 
This gives m = 300 g cm~ 2 s _1 ss 0.003mEdd, below the 
accretion rates of persistent LMXBs or most transient 
LMXBs in outburst. However, separation of light and 
heavy elements is something that should be considered 
during quiescent periods in transient accretors (Brown 
et al. 2002) or at low accretion rates (Peng, Brown, & 
Truran 2007). 
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2.3. Compositional buoyancy and convection 

After the solid particles form and sediment out, the 
fluid left behind is lighter than the fluid immediately 
above it and so will have a tendency to buoyantly rise. 
This is counteracted by the thermal profile, which is sta- 
bly stratified in the absence of a composition gradient 
such that a rising fluid element will be colder than its sur- 
roundings and will tend to sink back down. A measure of 
the buoyancy is the convective discriminant A, which is 
related to the Brunt- Vaisala, frequency A 2 = — gA (Cox 
1980). For a two-component mixture in the ocean, we 
can write (Bildstcn & Cumming 1998) 
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here, A is the mass fraction of the lighter element (Fe 
or O in the examples above), xx = <9hiP/<91n A|t, p , 
X P = d\nP/d\np\ X: T, XT = d\nP/d\nT\ x , P , and 
the temperature and composition gradients are V = 
-H P {d\aT/dr) and Vx = -H P (d\nX/dr). The adi- 
abatic gradient is taken at constant entropy S and com- 
position: V ad = -Hp(dln T/dr\ s ^x)- Note that xx, XT, 
and Xp are au positive quantities. If A < or TV 2 > the 
ocean is stable to convection. For example, if the compo- 
sition is uniform so that Vx = 0, stability to convection 
requires the familiar condition V < V a d- Only one term 
describing the variation of the composition is needed in 
Eq. (J7J since we consider a two-component mixture. We 
generalize to more than two species in the next section. 

In steady state the composition profile is lighter with 
increasing depth: Vx > 0. Such a profile will not lead to 
convection as long as the gradient is small enough that 
the destabilizing effect of the composition profile is com- 
pensated by the thermal buoyancy represented by the 
first term in Eq. ([7|. The maximum stable composition 
gradient is 
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where we assume that the large thermal conductivity in 
the ocean due to the degenerate electrons results in an 
almost isothermal profile V <C V a d (see SJ). 

As accretion continues, light elements are continually 
deposited at the base of the ocean and must be trans- 
ported upwards by convection. We expect therefore that 
the composition gradient will adjust to be close to but 
slight greater than Vx,max so as to result in the required 
convective flux of composition Fx — th{Xo — At), where 
Ao is the incoming composition and Xb is the composi- 
tion at the base of the ocean (that results in freezing of 
solid with mass fraction Ao). 

We can estimate Vx — Vx,max in steady state and the 
corresponding convective velocity using mixing length 
theory. The acceleration of a fluid element is gA, giv- 
ing a convective velocity 
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where l m is the mixing length and we again assume 
V = for simplicity. After moving a distance l m , the 
mass fraction differs from its surroundings by an amount 
(l m /Hp)XVx, implying that there is a flux of composi- 
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Setting this equal to the steady state flux at the base of 
the ocean Fx — m(X — A&) = pu accr (A — A&) gives the 
maximum convective velocity 
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In the deep ocean where pressure is dominated by de- 
generate electrons, all factors on the right hand side of 
Eq. ()llj) are of order unity except for xx ~ XpO- ~ 
Y Se /Y e ) ~ 0.1, where Y Sc = 34/79; and X t « 
(T/Pe)(dPi/dT) ~ 10k B T/[(Z)E F ] (e.g., Hansen & 
Kawaler 1994), where P e is the electron pressure, Pi 
the ion pressure, and Ep is the Fermi energy exclud- 
ing rest mass. For an ocean temperature of 3 x 10 s K 

and Ep — 5.1p^ 3 Ye^ 3 MeV (the electrons are relativis- 
ts), xt ~ 0.001. Therefore, w C onv is about two orders of 
magnitude larger than the accretion velocity. Given that 
chemical separation is being driven by accretion, it may 
be surprising that v conv is much larger than w accr . The 
reason is that only a very shallow composition gradient 
can be tolerated in the deep ocean where kpT <C Ep, 
requiring t> conv ^> v accr to transport the required flux of 
composition. 

Comparing Eqs. ([9]) and (fTTj) . we see that Vx — 

Vxmax < (^accr/c,) 2 [xpXx(A -A fc ) 2 /V 2 dX |A 2 )] « 

1, where c s — {gHp) 1 / 2 is the sound speed. In other 
words the convective velocities needed to transport the 
flux of light elements through the ocean are very sub- 
sonic, implying the composition gradient is extremely 
close to the marginally stable gradient Vx ~ Vx.max- 
This is analogous to efficient convective heat transport 
for which V ~ V a d- 

Another piece of physics that could potentially play a 
role is diffusion. Inwards diffusion of composition will oc- 
cur down the composition gradient Vx, and in principle 
if efficient enough could mediate the need for convection. 
The diffusive flux is much smaller than the convective 
flux, however. To see this, we use Eq. ([3]) for a selenium 
OCP: 
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The diffusion time across a composition scale height 
l x = Hp/Vx, where the pressure scale height Hp — 

1900pg /3 (re/0.43) 4 / 3 ( 5 i 4 /2.45)- 1 cm is 
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This is much longer than the accretion timescale given 
by Eq. ([2"]). or the convective turnover time given by 
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implying that microscopic diffusion does not play a signif- 
icant role in transporting light elements across the ocean. 

We also note that our conclusion that the ocean is con- 
vectively unstable depends on the slope of the liquidus 
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curve in the phase diagram. In all of the cases we have 
considered, the composition profile that would exist if 
the ocean followed the liquidus curve as fluid elements 
move to higher pressure is too steep to be maintained by 
thermal buoyancy. If this were not the different 
picture would result: the composition in the ocean at a 
given pressure would correspond to the composition of 
the liquidus curve at each depth, and a steady hail of 
solid particles would fall through the ocean to the solid 
crust at the base. Given the phase diagrams for two- 
and three-component plasmas we have calculated, how- 
ever, this situation appears not to arise. 

We therefore expect that the ocean adjusts its com- 
position from Xq at the top to at the base with a 
mixing zone in which the gradient is Vx ~ Vx,max- Be- 
cause this gradient is very shallow in the deep degenerate 
part of the ocean (Vx,max ~ 10kBT/[(Z)Ep]), we ex- 
pect the mixing zone to have substantial thickness. We 
investigate this in the next section with detailed mod- 
els of the steady-state ocean, and find that the entire 
ocean is expected to undergo mixing all the way up to 
the layer of light elements that supplies the ocean with 
new material through nuclear burning. The shallow gra- 
dient also means that there is a nearly uniform compo- 
sition throughout the bulk of the ocean, which is set by 
the phase diagram at the freezing point. 



3. ILLUSTRATIVE MODELS OF A CONVECTIVE 
OCEAN 

In <j2]we argued that crystallization and sedimentation 
would drive a convective instability at the base of the 
ocean. We now calculate detailed but illustrative models 
of the steady-state ocean based on this picture. We use 
mixing length theory to calculate the convective velocity, 
keeping the mixing length as a free parameter; whether 
mixing length theory is appropriate for compositionally- 
driven convection is an open question. As we discuss 
further in $51 we also neglect other hydrodynamical cir- 
culations (e.g. Eddington- Sweet circulation driven by ro- 
tation), the effects of rotation or magnetic fields on the 
convection itself, or other transport mechanisms such as 
turbulent mixing that could lead to transport of compo- 
sition through the ocean. 

We first allow the ocean to extend to arbitrarily low 
density; i.e., we first neglect the light element layer which 
overlies the ocean ( ^3.11) and then show how the light 
element layer can be incorporated self-consistently into 
a steady-state model of the ocean ( §3.2p . The models in 
£ 33.11 and ^3.21 are isothermal; in S|4]we consider the effect 
of the mixing on the thermal profile. 



3.1. A first model of the ocean 

We first calculate the steady state in the ocean for 
the case where nuclear reactions are unimportant in the 
convection zone. We show results for the specific two- 
component mixtures discussed in $2j but for generality 
we keep the number of species arbitrary in the deriva- 
tions that follow. The continuity equation for the flow of 
species i is 

+ V ^ccr ' VI, = - -V ■ (pVco^DXi) , (15) 

at p 
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for an inward-convecting blob, and the convective veloc- 
ity from mixing length theory is given by (e.g., Kippen- 
hahn & Weigert 1994) 



gH P 



XT(V - V ad ) 



n-l 

£ 

i=l 



(17) 



where n is the total number of chemical species in the 
ocean, Xi is the mass fraction of species i, and the pa- 
rameter £ = l m /Hp is the ratio between the convection 
mixing length l rn and the scale height. In mixing length 
theory the value of £ is highly uncertain; we assume that 
£ = 1 in this paper, but £ could be an order of magnitude 
or two smaller than this. 

In steady state and for vertical flow in plane-parallel 
geometry (a good assumption since the ocean is thin com- 
pared to the stellar radius), we have 

dXi £ d . . . . 

faccr-l— = 7T~T (P^conv A, V X t ) • (18) 

dr Ip dr 



Integrating from the top of the ocean to some depth 
within the ocean, using the fact that pw a ccr = in is a 
constant, we find 

Waccr(Aj - Xifi) = V conv ^X 2 V Xi ■ (19) 

Note that at the top of the convection zone, when v CO nv = 
0, we have X^ = X i for each Xi. Rewriting Eq. (|19|) 
gives 
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which can be integrated for each species i — 1, 2, . . . , n— 1 
(the composition of species n follows from the constraint 

The n — 1 differential equations (|20j) are coupled be- 
cause the convective velocity depends on a sum over all 
species. To obtain an expression for the convective ve- 
locity, we multiply Eq. (|19j) by xx, and sum over species 
to obtain 
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We argued in S}2] that small convective velocities v conv <C 
c s are required to transport the composition flux. There- 
fore, 

n-l 

Exx.Vx* «XT(V ad -V) (22) 
i=i 

[cf. Eq. ©], and we can replace the left-hand side of 
Eq. (|2T|) with x-r(V a d — V). Rearranging, we find an 
expression for v conv : 



1 



£xt(V 



ad 



5^ Xx, 



Xi — X, 



i.O 



i=l 



Xi 



(23) 



As a check we see that for two species (n = 1), Eqs. (|20p 
and (23) give din X/ din P = (xr/xx)(V ad - V) which 
is the marginally stable composition gradient ( ^2.3|) . 




Fig. 2. — (Color online) The composition profile in the convection 
zone of a T = 3 X 10 s K ocean composed of 56 Fe, 79 Se (solid line) 
and 16 O, 79 Se (dashed line). The hydrogen/helium burning layer 
is neglected. 



Fig. 3. — (Color online) The accretion velocity for m = 3 X 
10 4 g cm -2 s _1 , and the convective velocity within the convection 
zone of a T = 3 X 10 s K ocean composed of 56 Fe, 79 Se (solid line) 
and 16 Q, 79 Se (dashed line). 



To obtain a solution, the location and composition of 
the base of the ocean are specified, and then Eq. (|20[) is 
integrated outward for each JQ . The top of the convec- 
tion zone is located at the point where the composition 
matches the incoming composition, X, t — X i0 . In Fig. [5] 
we present our results for the 56 Fe- 79 Se and 16 0- 79 Se sys- 
tems described in <21 The ocean temperature is taken to 
be T = 3 x 10 8 K. To calculate the various V's and 
X's, we assume that the electron pressure is given by the 
fitting formula of Paczyhski (1983), and following Paper 
I we include Coulomb corrections for the ions using the 
free energy from DeWitt & Slattery (2003). Note that 
the equations used in this section to calculate the compo- 
sition profile are all independent of £ ; only the convective 
velocity depends on this value. Therefore, the potentially 
large error introduced by our choice of £ does not affect 
our results. 

The solutions show the expected behavior, that at large 
depths where fc^T -C Ep, the composition gradient is 
very shallow, but as the integration continues further 
upwards, the gradient steepens as Ep drops. In each 
case, the ocean is substantially enriched in light elements 
throughout most of its mass. Figure |3] compares the con- 
vective velocity with the accretion velocity for these two 
cases, assuming £ = 1 and = 3 (i.e., an accretion rate 
of 0.3mEdd)- Towards the top of the ocean, where the 
composition gradient becomes significant (Ep becomes 
comparable to fc^T), the convective velocity drops to- 
wards the accretion velocity. 

3.2. Including nuclear burning of light elements 

Figure [5] shows that in the oxygen- selenium two-species 
model at T = 3 x 10 s K, the convection zone extends into 
the hydrogen/helium burning layer (p ~ 10 5 -10 7 g/cm 3 ). 
This effect becomes more pronounced at lower tempera- 
tures, as might be expected when m < 0.3rh,Edd: empir- 
ically the density at the top of the convection zone goes 
as T 2 - 3 for typical ocean temperatures, in both the O- 
Se and the Fe-Se models. Therefore, an accurate model 
of the entire convection zone must include the effects 
of nuclear reactions. As a first approximation, we use 



the following simplified model for the burning layer: At 
p = 10 g/cm 3 the neutron star is composed only of hy- 
drogen and helium, with Ah — 0.7 and Xn e = 0.3. The 
relative abundance of hydrogen and helium is maintained 
as these elements undergo nuclear reactions and convec- 
tive mixing, such that An/Ane = 0.7/0.3 is a constant 
throughout the ocean. This mixture burns at the triple- 
alpha rate (Hansen & Kawaler 1994) 



R: 



-3a 



e -44/T s 



,-1 



(24) 



where Xi is the mass fraction of the H/He composite 
with Ah = 0.7Ai and An e = 0.3Ai. Each gram of the 
mixture burns to Xf UTn grams of each chemical species 

i > 1 (so that X ? Um = !)■ 

We use this model for the burning layer to calculate 
the steady state in the ocean for the case where nuclear 
reactions are important in the convection zone. The con- 
tinuity equation for each species becomes 

^+V aC cr'VIi = --V-ipVcoavDXil-XiRaa (25) 

at p 
for i — 1 and 

^ + v accr • VI, = - iv • (pv conv Z?A,) + X^ n X,R 



dt 



P 



3o 



for i = 2, 3, . . . , n. In steady state we have 
dX 1 £ d 



(26) 



dr 2p dr 



(pV co nvXiVxi) + X 1 R 3a (27) 



and 



«accr^ = (PWA,V X J - Af^A^a . (28) 

dr lp ar 



Above the convection zone (where v conv = 0) we use 
dX l 



dlnP 



-XiR 3a Hp 



(29) 
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to solve for X\ as a function of P, and then 



d\nP 



Tyburn y E> TJ 



or 



X = (i - x^x^™ 



(30) 
(31) 



to find the Xj, i > 1 values. These equations are solved 
inward from p = 10 4 g/cm 3 , where we set X\ — 1 and 
Xi = for i > 1. Within the convection zone we use 
the following method to find the X^'s: We first spec- 
ify the location of the top of the convection zone, Pq. 
The Xifl values used in the equations below are then 
found by solving Eqs. (f2"Tj|) and (f3~Tj) down to this point. 
The value of Pq is not known a priori, so the Xj's are 
found for a given Pq and then the value is varied until 
the composition changes smoothly across the convection 
zone boundary. We define a new variable, W, such that 



dW 
dr 



X X R, 



3 a 



(32) 



and W — at the top of the convection zone (when Xi - 
Xifi). Then the equation for the convective velocity is 



«-i y v 

Xx, ~ - W 



i=l 



X; 



n-1 N 
V"burn 



^conv £ 



^Xx.V^, (33) 



»=i 



or 



Xj — X; o 




The n coupled differential equations to solve are 

dX\ Waccr 



dlnP 



1,0 



dlnP 



(X, - X 4 , - WX 



W): 



burn\ 



(34) 



(35) 



(36) 



for i = 2, 3, ...,n - 1; and Eq. (|3"5]). We solve these 
equations through iteration, first guessing the Xi values 
at every point in the convection zone, then using these 
values to find W, and so on, until convergence is reached. 

In Fig. [4] we present our results for the 56 Fe- 79 Se 
and 16 0- 79 Se systems discussed in Sj2j again we choose 
7714 = 3, £ = 1, and an ocean temperature T = 3 x 10 8 K. 
For comparison we also present the Xp e , Xq, and X$ e 
profiles when the hydrogen/helium burning layer is ne- 
glected f ^3.1[) . These profiles are nearly the same as 
the profiles when the burning layer is included, except 
at the top of the convection zone where the latter pro- 
files are very steep. The change in the slope of the pro- 
files from §3.11 to §3.21 can be understood from Eq. (|22|) : 
This equation states that within the convection zone 

J27=i Xx ( Vx< follow the adiabatic gradient, but places 
no such restriction on the individual X profiles. There- 
fore, at the top of the convection zone, the Fe, O, and Se 
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Fig. 4. — (Color online) The composition profile within and di- 
rectly above the convection zone ofaT = 3xl0 8 K ocean composed 
of 56 Fe, 79 Se (top panel) and 16 O, 79 Se (bottom panel). A layer 
of hydrogen and helium is placed on top of the ocean and nuclear 
reactions are crudely modeled. For comparison the Xp e , Xq, and 
Xg e profiles when the burning layer is neglected (Fig. , and the 
. X"h + Xne profiles when convection is turned off [Fqs. H29H and 
(131 H i are also plotted, as thin lines with the same patterns as the 
profiles from the full calculation. The convection zone begins at 
the point where the two Xn + X^ e curves diverge, at a density of 
7.2 X 10 6 g/cm 3 for the F-Se system and 2.6 X 10 6 g/cm 3 for the 
O-Se system. 



composition profiles rise sharply with increasing density 
because of the correspondingly large drop in the H-Hc 
profile. The burning layer acts as a barrier to convec- 
tion, such that the convection zone ends abruptly at the 
base of this layer. 

Finally, in Fig. [4] we present Xh + Xh c profiles when 
convection is ignored [such that Eqs. (f29| and (|3"T|) ap- 
ply across the entire ocean]. The shape of these pro- 
files results from the balance between accretion driv- 
ing the H-He mixture deeper and nuclear reactions con- 
verting H-He into heavier elements as it travels inward; 
the profiles can be estimated using Eqs. ([2]) and ([24)) : 
1— i a ccr-R3a — Xi ^> Xi cx p 5 / 3 for Xi < 0.1. In addition 
to the sharp drop at the top of the convection zone, the 
profiles when convection is included differ in that they fall 
more slowly with increasing density, due to their larger 
inward velocity; since w CO nv oc p~ 2 / 3 in the deep ocean, 
Xi oc p 3 / 2 for Xi < 10~ 3 . The X H + X Hc profiles there- 



fore cross at some depth in the ocean; i.e., convection 
causes some hydrogen-helium to mix into deeper layers. 
While this mixing will change the amount of burning 
at each depth in the convection zone, we do not expect 
the thermal profile in this region to change, because the 
amount of helium in the convection region is very small 
(less than 10% by mass of the total helium in the star, 
and less than 1% of the total mass at any depth in the 
convection zone; see Fig. [4j. 

4. EFFECT OF MIXING ON THE THERMAL 
PROFILE 

So far we have taken the ocean to be isothermal, which 
is a good first approximation for the bulk of the ocean 
where thermal conductivity efficiently transports heat 
(Bildstcn & Cutler 1995). In this section, we discuss 
the effect of convective mixing on the thermal profile 
and evaluate the heating associated with the transport 
of light elements upwards through the ocean and the 
steady-state thermal profile. 

When mixing is unimportant, the main contribution to 
the heat flux in the ocean comes from energy released by 
nuclear reactions in the crust. These reactions release an 
energy of « 1.7 MeV per nucleon (e.g., Haensel & Zdunik 
2008), of which Qb ~ 0.1-1 MeV per nucleon flows out- 
ward through the ocean, depending on the accretion rate 
(Brown 2000). This component of the heat flux can be 
written as 

Fcrust = Qb — ■ (37) 
m p 

The temperature gradient required to carry the heat flux 
is given by 

FH P , . 

V = , (38) 



where 



K = 



KT 



Aac 

3pKrad 



K, 



cri 



(39) 



is the total thermal conductivity including both radia- 
tion and conduction contributions, a is the radiation con- 
stant, c is the speed of light, K ra d is the radiative opacity 
(e.g., Schatz et al. 1999), K c< \ = n e k 2 B T '/ '(m*J/ c ) is the 
thermal conductivity, = m e 4- Ep/c 2 is the effective 
mass of the electrons, and v c is the collision frequency. 
In the ocean, the relevant collision frequency is that be- 
tween electrons and ions: v c — 4e 4 m*A/ (3irh 3 )(Z 2 ) / (Z) , 
where A is the Coulomb logarithm (Yakovlev & Urpin 
1980; Schatz et al. 1999). For densities p > 10 7 g/cm 3 
the electrons are degenerate and relativistic, and the con- 
tribution of radiation to K is negligible; in this case we 
find 



V ~ 0.03 m 4 



Qi, 



0.1 MeV J V 3 

2 



x 



mm)fn±y\ (40) 

30 J V2.45/ ' V ' 



where in the above we have set A = 1. We see therefore 
that over most of the ocean a shallow temperature gradi- 
ent is sufficient to conduct the heat flux from the crust, 
and in particular V <C V a( j ~ 0.35 as we assumed in fj3j 
Note that this is not true for p < 10 7 g/cm 3 , where the 



electrons are nonrelativistic. In that regime Hp/K c & oc 
p -1 / 3 , such that at low densities FH P /(K cd T) > V ac i. 
The radiation contribution to the total K becomes im- 
portant at these densities, partially offsetting the effect of 
a low K c d- Nevertheless, under certain conditions (e.g., 
rh > mEdd and T < 10 9 K) V > V a d at the top of 
the ocean. It is unclear what happens in that case; a 
complete understanding may require a time-dependent 
calculation (see fj5|. 

Thermal conduction also readily conducts the latent 
heat away from the ocean floor. We can estimate the 
latent heat release using the equations from Paper I for 
the free energies of the liquid and solid states [e.g., equa- 
tions (23) and (24) of that paper]; the latent heat released 
on freezing is of order k B T/(A) « 0.1 T 8 [60/(A)] keV 
per nucleon. Much smaller than Qb, the latent heat is re- 
moved with only a small temperature gradient. This con- 
trasts with, e.g., freezing of solid material in the Earth's 
core, in which the temperature gradient required to con- 
duct the latent heat away is V > V a d, so that the latent 
heat release drives thermal convection (Stevenson 1981). 

When chemical separation drives convection in the 
ocean, we must also consider the convective heat flux 
in addition to the conductive heat flux. In mixing length 
theory the heat flux is 

Fconv = |pv conv c P T(V - V ad ) , (41) 

where cp is the heat capacity. Using Eq. (|2"5|) for the 
convective velocity together with the fact that pVn Ccr = 
rh, we find 

Fconv — E conv T (^^) 



with 



F r , 



cpTrh ^ X t - Xt.o 



XT 



Xi 



(43) 



Here r is the radial unit vector. Note that F Gonv > 
so that the convection transports heat inwards. This is 
because, unlike thermally-driven convection, V < V a a so 
that a fluid element displaced adiabatically outwards is 
cooler than its surroundings at its new location. 

We show the convective flux as F conv m p /rh with rh = 
3 x 10 4 g cm~ 2 s _1 in Fig. [5l In the ocean, the heat 
capacity is set by the ions; using the internal energy ex- 
pansion from DeWitt & Slattery (2003) we find cp ~ (1- 
A)k B /[(A)m p ] giving c P T/m p ~ (0.1-0.6)Ts[60/(A)] keV 
per nucleon, where the range of values is across the 
depth of the ocean. Although cpT/m p is much smaller 
than Qb, the convective flux has an additional factor of 
Xx/xt Ep /(kpT). This gives an extra factor of 10- 
100, so that the final convective flux is - 10~ 2 -10 _1 MeV 
per nucleon, which can be comparable to Qb- For the O- 
Se ocean in Fig. [5J the convective flux is ps 0.2 MeV per 
nucleon at the base of the convection zone. For Fe-Se, the 
smaller contrast between the heavy and light elements 
gives a much smaller flux. For two species, the compo- 
sition flux is Fx = rh(X — Xq), so that we can write 
F com = F x c p Txx/(Xxt) = F x c P {d\nT/d\nX\ P , p ). 
The convective heat flux corresponds to the rate of 
change of internal energy given the flux of composition 
at each depth [compare equation ([2]) of Montgomery et 
al. 1999]. 




Fig. 5. — (Color online) The convective flux for a fn 



10 g cm 



(solid line) and ia O, 79 Se (dashed line). 

Figure [3] shows that the inwards convective flux in- 
creases with depth in the ocean dF conv /dy > 0. This 
implies that there is a local cooling at each depth in the 
ocean, and so while it is a good approximation to assume 
the ocean is isothermal when calculating the convective 
velocity and composition profile of the mixed zone, in 
fact the isothermal models we presented in SJ3] are not in 
a thermal steady-state. To see what the thermal steady- 
state must look like, we write down the entropy equation 



3 x 

3 X 10 s K ocean composed of 5b Fe, 79 Se 



T D 1 = dF 
Dt dy 



(44) 



where s is the specific entropy, e is the sum of all sources 
and sinks of heat (latent heat release, nuclear reactions, 
neutrino cooling, etc.), and D/Dt = d/dt + fnd/dy is 
the total derivative for a fluid element. The heat flux F 
is the sum of convective and conductive heat fluxes. In 
steady state (ds/dt = 0) the left hand side of Eq. (|4"4")l is 
given by (e.g., Brown & Bildsten 1998) 



T Ds = c_pTm 
Dt y I 



^ n— 1 



(45) 



much smaller than the dF/dy terms on the right hand 
side since in the convective region Y17=i Xx^ Xi ~ 
XT(V a d — V) [Eq. (|22]) ] (the convection zone is very close 
to adiabatic). Neglecting any contributions to e in the 
ocean, e.g. nuclear reactions or the small latent heat 
release at the boundary, we must have dF/dy = in 
steady state, or an outwards conductive flux -F C ond = 
-Fcrust + -Fconv The picture we arrive at is that start- 
ing with an isothermal ocean, the inwards convection of 
heat steepens the temperature gradient until the out- 
wards conductive flux balances the convective flux. 

We see therefore that the effect of the mixing is to add 
a contribution to the conductive flux in the deep ocean 
of as much as 0.2 MeV per nucleon for an O-Se mixture. 
At high accretion rates fn > 0.1 WEdd, this is comparable 
to or larger than the flux from the crust Qb ~ 0.1 MeV 
per nucleon (Brown 2000; at lower accretion rates fn ~ 
0.01 mEdd, a larger fraction of the heat released in the 
crust flows outwards, and Qb ~ 1 MeV). We show in 



Fig. 6. — (Color online) The thermal profile for a fn = 3 X 
10 4 g cm -2 s _1 ocean composed of ls O, 79 Se and with T = 
2 x 10 s K at p = 8 x 10 5 e/cm , both when the convective flux is 
included in the total heat flux F (solid line), and when it is ignored 
(i.e., the total heat flux is due only to the outward flux from the 
crust; dashed line). 

Fig.[6]the thermal profile with and without the convective 
heating included, for an O-Se mixture with 7714 = 3 and 
T = 2 x 10 8 K at p = 8 x 10 5 g/cm 3 . 

5. DISCUSSION 

We have explored the consequences of chemical sepa- 
ration in the ocean of accreting neutron stars. Given the 
rapid timescales for nucleation, growth, and sedimenta- 
tion of solid particles, fluid elements that are lighter than 
their surrounding are continually being released at the 
base of the neutron star ocean. Using a mixing length 
model for convection, we modeled the resulting mix- 
ing zone. The conclusions are that the entire ocean is 
mixed, with a composition gradient Vx ~ XT^ad/xx ~ 
100k B T/[(Z)E F ] that is very shallow in the bulk of the 
ocean where the electrons are relativistically degenerate 
(Ep 3> ksT). The composition profile is therefore rather 
uniform and enhanced in light elements compared to the 
composition produced by nuclear burning of the accreted 
light elements at the top of the ocean. For example, for 
the largest charge ratio we considered, a mixture of O 
and Se, only 2% oxygen by mass was added to the ocean 
by nuclear burning, but the ocean itself is enriched in 
steady state to almost 40% by mass of oxygen (Fig. (2), 
set by the phase diagram for the O-Se mixture. 

Accompanying the outwards transport of light ele- 
ments is an inwards transport of heat. In £|4l we argued 
that the ocean will evolve to a steady-state in which the 
convective heat flux is balanced by an outwards conduc- 
tive flux. The effect of chemical separation is therefore to 
add a contribution to the outwards conductive heat flux 
in the ocean, steepening the temperature gradient. For 
the O-Se case, we find a heat flux of 0.2 MeV per nucleon 
at the base of the ocean, comparable to the flux expected 
from the crust. The heating is smaller for an Fe-Se ocean, 
which has a smaller contrast in composition between the 
two species. However, only a small amount of oxygen in 
the material entering the ocean is required to generate 
significant heating. In the O-Se case we considered, only 
2% of the incoming mass is oxygen, compared to 40% 
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following enrichment. Therefore, in an ocean consisting 
of many species, even a small amount of light element en- 
tering the ocean could lead to significant heating through 
its enrichment. We are currently investigating the phase 
diagrams of multicomponent mixtures to assess the level 
of heating in that case. 

An immediate application of our results is to super- 
bursts. Following Horowitz et al. (2007) and Paper I, 
we chose oxygen as the light element for our example 
in this paper, but we have calculated models with car- 
bon as the light element with similar results. Carbon 
enrichment through chemical separation at the base of 
the ocean, and subsequent upward mixing, could bring 
the carbon mass fraction Xc at the superburst ignition 
depth to within the required range Xc ~ 0.1-0.2 inferred 
from observations. Similarly, the heating associated with 
chemical separation will help to bring the temperature 
of the ocean up to the required ignition temperatures 
of 5-6 x 10 8 K. Questions still remain, however, such 
as whether this picture could match the observed recur- 
rence times. The hotter ocean could also drive a heat flux 
into the outer crust, explaining the inverted temperature 
gradient found by Brown & Cumming (2009) in cooling 
transients. Time-dependent calculations of the evolution 
of the ocean are in progress to address these issues. 

Compositionally-driven convection is important in 
cooling white dwarfs as well. The main difference be- 
tween convection in the cores of cooling white dwarfs 
and in the oceans of accreting neutron stars is that in 
the white dwarf case a given mass element crosses the 
liquidus curve and crystallizes because the star cools, not 
because the element is pushed inwards. The similarities 
between these two systems are numerous. For example, 
the convectively stable composition gradient in the white 
dwarf interior is Vx,wd ~ 0.03, nearly the same as the 
gradient in the neutron star ocean; the white dwarf crys- 
tallization rate, towd ~ 0.1 g cm~ 2 s _1 , is also similar 
(Mochkovitch 1983). Therefore, we expect that much of 
our analysis and calculation from the present paper is 
applicable to the white dwarf case. 

Another situation in which compositionally-driven 
convection is important is in the Earth's core, in which 
chemical separation at the boundary between the liq- 
uid outer core and solid inner core drives convection in 
the outer core and provides an energy source for the 
Earth's dynamo (Stevenson 1981). There are interesting 
differences between the Earth's core and the neutron star 
ocean, however. First, the thermal gradient V required 
to carry away the latent heat from the Earth's inner- 
outer core boundary is significant, and exceeds V a d, so 
that the latent heat drives thermal convection. Second, 
in the Earth's core V a d < Vl, where Vx is the liq- 
uidus temperature gradient (how the melting tempera- 
ture varies with pressure). This allows the boundary 
between the inner and outer core to be at the melting 
temperature, but the outer core with an adiabatic pro- 
file remains above the melting temperature. Therefore a 
fluid thermally convecting outer core is possible (Loper 
& Roberts 1977; Stevenson 1980; the opposite was pro- 
posed by Higgins & Kennedy 1971 in what became known 
as the "core paradox" ) . In the neutron star ocean, how- 
ever, V a d > Vl: Equation ^ shows that for fixed 
P oc T 4 so that Vl = 0.25; whereas using the results 



of DeWitt & Slattery (2003) to calculate Coulomb cor- 
rections we find V a d ~ 0.36 with a weak Ti dependence 
near the melting point. This means that the entire ocean 
could never be thermally convective and remain liquid. 

The steady-state convective velocity is larger than the 
accretion velocity by a factor (X^ — Xq) / (AbV x), where 
Xf, — Xq is the required change in composition across 
the ocean. Because the marginally-stable composition 
gradient in the ocean is very small, this factor is much 
greater than unity; in the deep ocean, this factor is of 
order 100. However, the corresponding velocity is still 
extremely slow such that additional processes that we 
have not considered here could change the picture we 
have put forward. For example, an extremely weak mag- 
netic field (only about 100 G at the base of the ocean) 
would have an energy density comparable to the kinetic 
energy of the flow. The dipole magnetic fields of ac- 
creting neutron stars in low mass X-ray binaries are be- 
lieved to be ~ 10 9 G and stronger horizontal field com- 
ponents would not be surprising, e.g. due to screen- 
ing currents (Cumming, Zweibel, & Bildsten 2001) or 
winding due to differential rotation and subsequent in- 
stabilities (Piro & Bildsten 2007). On a rotating star, 
Eddington- Sweet circulation in the ocean would have a 
timescale ~ ithcrm(^s/^) 2 where Q,b is the break-up 
spin frequency. This timescale could be hundreds of 
days for a thermal time of days at the base of the ocean 
and tt/Qs ~ 0.1, comparable to the convective turnover 
time. Further work is needed to investigate how rota- 
tional circulation or magnetic stresses would affect the 
models presented here. 

The steady state we have been discussing is one in 
which the unstable composition gradient in the ocean 
is balanced by the stable temperature gradient, such 
that the ocean is only marginally unstable to convec- 
tion. We find, however, that under certain conditions 
(e.g., m > rh-Edd) the top of the ocean is thermally unsta- 
ble to convection: the temperature gradient must satisfy 
V > Vad m order to conduct away the large flux from the 
crust. In this case, either the convective velocity must 
be very large at the top of the ocean, or the unstable 
temperature gradient must be balanced by a stable com- 
position gradient. The latter scenario is inconsistent with 
the two-component models of £)3.1l where the fraction of 
the light element increases with depth; but could be con- 
sistent with the multicomponent models of §3.21 if the 
fraction of hydrogen-helium drops faster than oxygen or 
iron rises. We are currently studying the time-dependent 
case to understand what happens when V > V a d- 

In addition to oceans with thermally-driven convec- 
tion, there are several other scenarios that are likely 
to lead to interesting time dependence of the mixing 
zone. For example, in transients in quiescence (see, e.g., 
Shternin et al. 2007; Brown & Cumming 2009) the mix- 
ing zone relaxes on a timescale comparable to the con- 
vective turnover time [Eq. (|14l) ]. about a month in the 
deep ocean. Mixing also occurs as the star cools after an 
accretion episode; the bulk of the ocean will solidify after 
a few days of cooling and will undergo chemical separa- 
tion and further convection. Either of these effects could 
lead to late-time energy release, and provide a further 
observational test of the conditions in the ocean. 
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